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ABSTRACT 


An  efficient  method  of  computing  structural  response  of  multi-story  nonlinear 
base  isolated  buildings  for  a  given  seismic  event  is  presented.  Using  a  recursive  block- 
by-block  integral  equation  formulation  (RBBIEF)  solution  to  the  governing  nonlinear 
Volterra  integral  equation,  structural  base  motion  acting  through  an  arbitrary  number  of 
nonlinear  base  isolators  can  be  computed  quickly  and  accurately.  The  general  solution  to 
the  governing  nonlinear  Volterra  integral  is  formulated  and  subsequently  converted  into 
code  using  MATLAB.  The  procedure  incorporates  modal  properties,  computed  from 
conventional  finite  element  (FE)  techniques,  and  the  generated  MATLAB  programs  to 
solve  a  varying  set  of  multi-degree  of  freedom  structures  coupled  to  both  linear  and 
nonlinear  isolators.  Ultimately,  an  analysis  is  conducted  on  a  30-story  building  that  was 
overly  designed  using  the  1994  Load  Resistance  Factor  Design  [Ref.  1]  and  the  1994 
Uniform  Building  Codes  [Ref.  2]  for  earthquake  loading.  The  method  demonstrates  that 
the  Volterra  integration  scheme  in  the  time  domain  is  very  effective  and  efficient. 
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I.  INTRODUCTION 


Seismic  base  motion  is  an  important  design  factor  in  constructing  and  retrofitting 
structures.  The  extensive  base  displacements  and  resultant  shear  forces  caused  by  large 
magnitude  earthquakes  are  major  contributors  to  the  eventual  failure  of  the  structure. 
Isolating  structures  from  seismic  events  is  currently  accomplished  through  the  use  of 
passive  isolators. 

Originally,  it  was  the  engineer’s  goal  to  reduce  the  structural  displacements 
through  isolation.  However,  current  techniques  place  equal  importance  on  increasing  the 
energy  dissipation  in  an  isolation  device.  An  increase  in  the  area  contained  within  the 
device’s  hysteresis  loop  is  directly  related  to  the  energy  dissipation  provided  by  the 
device.  The  plastic  properties  of  passive  hysteretic  type  isolators  are  mathematically 
modeled  as  a  complex  nonlinear  equation  [Ref.  3],  Modeling  and  analysis  using  finite 
element  (FE)  techniques  for  nonlinear  isolators  then  becomes  equally  complex  and 
computationally  demanding. 

Several  Commercial-Off-the-Shelf  (COTS)  programs  such  as,  DDRAC-3D,  3D- 
BASIS,  DRAIN-2DX,  SAP,  ANSR,  N-PAD,  LPM,  and  others  have  their  respective 
limitations  and  advantages  as  outlined  in  [Ref.  3].  The  primary  solution  methods  for 
nonlinear  base  isolation  of  structures  are  the  pseudo-force  iterative  method  or  the 
nonlinear  tangent  stiffness  method  [Ref.  3].  However,  the  incorporation  of  3-D  motion 
and  discrete  time  histories  in  conjunction  with  n-nonlinear  isolators  to  the  above- 
mentioned  FE  programs  can  be  computationally  exhaustive  and  or  not  practically 
possible. 
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Current  FE  programs  can,  in  seconds,  compute  normalized  mode  shapes  and 
natural  frequencies  of  several  thousand  degree  of  freedom  systems.  The  structure’s 
modal  properties  are  used  by  many  of  the  above  programs  in  order  to  compute  transient 
response  or  frequency  response.  Unlike  the  methods  of  solving  the  differential  equations 
of  motion  above,  we  describe  herein  a  method  that  solves  the  transient  response  directly 
using  the  Volterra  integral  form.  This  method,  referred  to  as  time  domain  synthesis,  is 
presented  in  [Ref.  4]  and  [Ref.  5],  and  is  a  numerically  accurate  and  fast  solution  for 
structural  response  due  to  base  motion  excitation  through  nonlinear  isolation. 

The  addition  of  nonlinear  memory-type  isolators  to  the  time  domain  synthesis  is 
presented  herein.  Currently,  the  synthesis  is  solved  using  the  (RBBIEF)  method.  The 
RBBIEF  method  obtains  its  significant  efficiency  from  two  characteristics  of  the 
formulation.  First,  a  stmcture  is  described  by  impulse  response  functions  calculated  only 
for  those  physical  DOF  of  interest.  Those  DOF  for  which  response  is  not  required  can  be 
omitted  with  no  loss  of  accuracy.  Second,  the  block-by-block  convolution  divides  a  large 
time  history  into  several  “blocks,”  and  the  governing  Volterra  integral  equation  is  solved 
in  a  block-by-block  manner  wherein  only  one  block  undergoes  iterative  convolution  at  a 
time.  The  converged  forces  of  synthesis  from  previously  iterated  blocks  are  used  in 
conjunction  with  the  results  of  the  current  iteration  to  calculate  the  response  for  the 
current  block.  To  incorporate  memory-type  isolators  in  the  synthesis,  the  isolator 
equations  are  solved  using  standard  numerical  solution  techniques  for  ordinary 
differential  equations.  This  solution  of  the  isolator’s  governing  differential  equation  is 
incorporated  into  the  RBBIEF  method.  This  revised  synthesis  method  provides  fast  and 


accurate  solutions.  Ultimately,  the  combination  of  quick  eigenvalue  solutions  of  large 
degree  of  freedom  structures,  modal  reduction,  and  a  block-by-block  recursive  solution 
of  the  governing  integral  equation  of  the  retained  nodes  that  are  seismically  excited,  can 
lead  to  significant  timesavings  and  accuracy. 
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II.  EARTHQUAKE  DESIGN  AND  NONLINEAR  ISOLATORS 


A.  STRUCTURAL  DESIGN 

Current  design  techniques  and  survivability  rely  upon  the  structure’s  total  strain 
energy  capacity  as  it  plastically  deforms  during  a  seismic  event.  This  “hard  working” 
process  is  accounted  for  by  the  Uniform  Building  Code,  UBC.  In  Figure  1  below,  the 
UBC’s  required  seismic  coefficient  is  significantly  underestimating  the  various 
earthquakes  presented. 


Figure  1.  UBC  Seismic  Acceleration  Graph  From  Ref.  [3] 

Additional  varying  normalized  response  spectra  plots  similar  to  Figure  1  are 
produced  in  [Ref.  3].  The  vertical  axis  outlines  the  seismic  coefficient  and  is  used  in 
various  UBC  equations,  the  horizontal  axis  represents  structural  periods.  The  apparent 
differences  between  the  UBC’s  required  protection  as  seen  in  Figure  1  and  an  arbitrary 
seismic  excitation  is  where  passive  base  isolation  plays  it’s  important  role.  In  order  to 
ensure  the  reduction  of  the  earthquake’s  overall  input  forces  using  passive  isolators,  an 
assumed  normalized  response  spectrum  is  chosen. 
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Existing  codes  outlined  in  Chapter  23  of  [Ref.  3]  require  that  the  existing 
structure  have  several  design  capacities  such  as  seismic  base  shear,  lateral  seismic  force, 
lateral  deflection,  overturning  moment,  and  torsional  moment.  As  part  of  the  seismic 
design  procedure,  the  structure  is  converted  into  modal  space  and  a  conservative 
participation  factor  is  chosen.  The  general  design  equation  is  presented  as: 


Replacing  \|/  with  V  for  shear,  A  for  drift,  8  for  deflection,  M  for  moment,  and  T  for 
torsion  into  Equation  (1)  above  to  arrive  at  the  respective  design  limits.  The  subscript 
(m)  represents  the  modal  values  and  (x)  the  story  level.  Details  of  modal  decomposition 
and  participation  can  be  found  in  numerous  vibration  textbooks  such  as  [Ref.  6]  and  will 
be  summarized  in  subsequent  chapters  herein. 

A  30  story  moment  resisting  frame  (MRF)  building  was  conservatively  designed 
using  the  Manual  of  Steel  Construction,  Load  Resistance  Factor  Design  (LRFD)  [Ref.l], 
the  Uniform  Building  Code  (UBC).  [Ref.2],  and  a  steel  structure  design  textbook  [Ref. 
7].  The  structure  has  an  inter-story  height  of  17  feet  and  a  total  building  height  of  510  ft. 
The  columns  were  spaced  25-ft  center  to  center  leading  to  an  overall  length  of  400  ft  at 
the  base.  A  total  of  256  columns  are  present  corresponding  to  256  base  isolators  at  the 
base.  Each  floor  was  designed  to  support  a  dead  load  of  150  psf  and  the  roof,  50  psf.  The 
beams  were  divided  as  follows:  the  roof  consisted  of  W36xl35  beams;  Floors  30-21 
W36xl50;  Floors  20-11  W36xl60;  Floors  10-2  W36xl70.  The  first  floor  was  designed 
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as  the  foundation  however,  when  the  structure  was  redesigned  for  isolation,  W36xl70 
beams  were  used  to  connect  the  columns  at  the  base  of  the  structure. 

Because  the  isolated  and  non-isolated  structures  are  two  distinct  eigenvalue 
problems,  due  to  the  free-free  nature  of  the  isolated  structure,  vertical  motion,  not 
previously  seen  in  a  clamped  structure,  is  no  longer  constrained.  The  addition  of  vertical 
displacements  at  the  base  is  known  as  uplift.  This  uplift  motion  can  increase  local 
stresses  as  well  as  decrease  the  designed  overturning  moment.  The  effects  of  uplift  can 
be  seen  at  some  or  all  frequencies  and  can  have  an  adverse  effect  on  many  types  of 
passive  isolators. 

B.  EARTHQUAKE  SAMPLING 

Seismic  events  occur  worldwide  and  on  a  frequent  basis.  Generally,  earthquakes 
who’s  magnitudes  are  greater  than  five  can  cause  significant  damage.  The  following 
uncorrected  ground  acceleration  plotted  in  Figure  2  is  indicative  of  earthquake 
magnitudes  of  6  and  above.  Various  time  history  data  can  be  obtained  through  several 
resources  such  as  [Ref.  8], 
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IMPERIAL  VALLEY  CA  1  07 D  10  IS  231  6 


Figure  2.  Imperial  Valley  Accelerometer  Graph  From  Ref.  [8] 

Without  formally  verifying  the  enormous  earthquake  time  history  data  for 
ergodicity,  many  of  the  earthquakes  observed  lasted  for  less  than  30  seconds  with  the 
majority  of  the  energy  focused  within  the  first  10-20  seconds.  Time  sampling  from 
various  monitoring  stations  is  generally  on  the  order  of  10'2  seconds.  Therefore  time  data 
histories  can  total  to  about  6000  data  points.  To  reduce  the  amount  of  data  and  speed  up 
the  subsequent  numerical  solutions,  correlations  can  be  calculated  similar  to  data  signal 
processing. 

The  equations  presented  in  this  section  are  outlined  in  Reference  [9]  and  are 
presented  here  in  summary.  This  direct  characterization  is  best  seen  in  the  spectral 
density  function  plot.  Beginning  with  the  auto-correlation  function.  Equation  (2),  of  the 
acceleration  time  history  a  set  of  characteristics  can  be  determined  from  the  earthquake 
time  history. 
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Rxx(r)  =  lim  ■—  J  x(t)x(t  +  T)dt 

(2) 

T  — 1  o 

R^xiO)  =  mean  square  value 

(3) 

Rxx(cc)  =  mean  value  squared 

(4) 
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Although  Equations  (3)  and  (4)  are  useful,  Equation  (2)  is  necessary  in  the 
computation  of  the  power  spectral  density  function.  Taking  the  Fourier  transform  of 
Equation  (2)  a  one-sided  spectral  density  function  is  developed,  Equation  (5).  This 
function  is  better  known  as  the  power  spectral  density  function  (PDF). 

G;tt(/)  =  2{^(T)e-2^T  (5) 

0 

Plotting  Equation  (5)  versus  frequency  will  yield  a  graphical  representation  of  the 
energy  distribution  of  the  system.  There  are  many  algorithms  that  take  advantage  of  time 
averaging  of  the  sample  data  to  primarily  suppress  bias  errors  in  the  estimate.  From  the 
PDF  an  appropriate  sampling  frequency  can  be  determined. 

To  avoid  statistical  sampling  errors  the  minimum  sampling  frequency  should  be 
equivalent  to  the  Nyquist  frequency. 

tnyq=2*  f resolution  (6) 

To  increase  the  resolution  and  to  avoid  aliasing  the  general  rule  of  thumb  for  the 
sampling  frequency  is: 

=10*  1nyq  (7) 

Given  the  value  of  Equation  (7)  and  the  data  obtained  from  plotting  Equation  (5) 
versus  frequency,  the  previously  mentioned  6000  data  points  can  be  dramatically 
reduced.  Using  MATLAB’s  power  spectral  density  function  on  the  El  Centro  earthquake, 
of  1940  time  history  data  reveals  it’s  distribution  plot  in  Figure  3.  From  Figure  3  it  is 
apparent  that  the  majority  of  the  energy  is  well  below  1  Hz.  Therefore,  using  a  sampling 
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frequency  of  10  Hz,  the  time  step  equivalent  of  .1  seconds,  reduces  the  proposed  6000 
data  points  by  a  factor  of  10  to  600. 


Figure  3.  El  Centro  1940  North  South  Power  Spectral  Density  Function  of  the 
Displacement  Time  History  Data 

C.  NONLINEAR  BASE  ISOLATION 

Since  many  spectral  density  plots  of  seismic  events  indicate  a  dispersed 
frequency  range,  viscous  dampers  aren’t  as  efficient  in  energy  dissipation  than  isolators 
that  are  frequency  independent.  Frequency  independent  isolators  maintain  their  unique 
hysteretic  loop  across  a  wide  range  of  frequencies.  The  general  topic  of  base  isolation  is 
well  researched  and  discussed  in  [Ref.  10]  and  is  broadly  presented  in  this  sub-section. 

The  general  nonlinear  equation  that  describes  general  frequency  independent 
isolators  is: 

F  =  -K0(  1  +  iS)(x  -  Mg)  (8) 
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Where  K0  is  the  dynamic  stiffness,  8  is  described  as  the  loss  factor,  x  the  base  motion. 


and  ug  is  the  ground  base  motion.  Graphically,  Equation  (8)  spans  the  hysteretic  loop  as 


shown  in  Figure  4. 


Restoring 

force 


Figure  4.  Hysteretic  Loop  of  an  Elastic-Pure  Plastic  Isolator 
The  process  by  which  the  hysteretic  loop  is  formed  is  through  cyclical 
deformation  along  a  material’s  stress  strain  curve.  Deformation  below  the  tensile  and 
compressive  yield  stresses  describes  typical  linear  elastic  springs.  Deformation  beyond 
the  initial  yield  stresses  results  in  plastic  deformation.  In  this  region  the  material  resists 
additional  stress,  due  to  strain  hardening,  leading  to  a  secondary  stiffness  rate.  In  order  to 
generate  the  hysteresis  loop  the  subsequent  unloading  from  the  plastic  region  occurs 
elastically  until  it  plastically  deforms  again. 

The  basic  groups  of  base  isolation  in  use  today  are  linear  base  isolators  and 
nonlinear  energy  dissipation  isolators.  Generally,  the  isolators  are  not  designed  to 
plastically  deform  but  to  achieve  the  similar  hysteresis  loop  through  hybrid  systems. 
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Typical  linear  base  isolators  are  elastic  springs  such  as  rubber  or  high  damping  rubber 
(HDR).  Typical  nonlinear  isolators  are  hysteretic,  visco-elastic,  friction,  or  material 
devices.  The  preferred  isolators  in  use  today  are  the  hysteretic  nonlinear  devices  and  the 
friction  devices. 

1.  Rubber  Isolators 

Rubber  used  by  it  self  is  categorized  as  a  linear  base  isolator.  Its  main  advantage 
as  an  isolator  is  it’s  dynamic  shear  modulus,  G,  remains  relatively  constant  over  a  wide 
range  of  frequencies  and  temperatures  [Ref.  11].  Additionally,  it’s  dynamic  damping 
coefficient,  5,  in  Equation  (8),  is  relatively  constant  at  frequencies  lower  than  50  Hz 
[Ref.  11].  Given  that  the  majority  of  the  seismic  energy  was  concentrated  well  below  50 
Hz,  the  use  of  rubber  isolators  is  desirable. 

A  nonlinear  isolator  is  produced  when  different  types  of  rubbers  are  used.  The 
resulting  combined  stress  strain  curve  for  the  composite  rubber  isolator  is  tri-linear 
[Ref.ll.p.  199].  At  moderate  shear  strains,  such  as  wind  loading  and  service  loading,  the 
effective  spring  stiffness  is  very  low.  At  larger  shear  strains,  seismic  loads,  the  stiffness 
can  increase  dramatically  until  it  decreases  to  the  third  stiffness  rate.  However,  since  the 
deformation  is  contained  within  the  elastic  regions  of  the  isolator  there  is  no  hysteretic 
loop  generated. 

Additional  materials  or  devices  built  into  a  rubber  isolator  can  produce  a 
hysteretic  effect.  The  plastic  deformation  rate  would  dominate  the  isolator’s  property  at 
larger  deformations  while  the  elastic  stiffness  of  the  rubber  would  control  the  subsequent 
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moderate  cyclical  loading  and  unloading  on  the  isolator.  The  effect,  during  seismic 


loading,  is  generally  modeled  as  a  bilinear  model  as  shown  in  Figure  5. 


Figure  5.  Bilinear  Hysteretic  Loop 

Where  Ki  and  K2  are  the  two  distinct  stiffness  constants  and  the  variable,  yp,  is 
the  yield  to  post  ratio  parameter. 

2.  Elastomeric  Isolators 

Elastomeric  bearings  are  engineered  from  synthetic  rubbers.  The  properties  of  the 
bilinear  hysteresis  loop  can  be  readily  designed  and  manufactured  with  the  use  of 
synthetic  rubbers.  These  devices,  known  as  High  Damping  Elastomeric  Bearings 
(HDEB),  are  frequently  used  in  both  structural  and  industrial  use.  A  comprehensive 
study  conducted  in  [Ref.  12]  validates  work  on  performance  and  efficiency  issues 
between  polymer  type  bearings  and  High  Damping  Rubber  Bearings  (HDRB).  The  only 
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disadvantage  of  polymer  synthetic  rubber  isolators  is  that  over  time  the  tensile  and  tear 
strengths  become  degraded. 

3.  Hysteretic  Isolators 

Not  including  the  HDEB  or  HDRB  devices,  hysteretic  bearings  are  developed 
using  metals  or  composites  such  as  the  lead  core  rubber  isolator.  By  varying  materials 
and  manufacturing  processes  a  wide  variety  of  hysteretic  curves  can  be  engineered. 
Modeling  hysteretic  systems  is  mathematically  challenging  due  to  the  nonlinear 
governing  equations.  An  accepted  mathematical  model  to  Figure  3  above,  is  the  Wen 
nonlinear  Equation  [Ref.  13].  The  equation  is  very  adaptive  to  various  hysteretic  type 
bearings.  Another  model,  not  used  in  this  paper,  is  known  as  the  Ozdemir  element 
[Ref.  14],  This  mathematical  model  is  based  on  similar  principles  as  the  Wen  element 
however;  it  incorporates  modifications  to  allow  for  Shape  Memory  Alloy  (SMA)  device 
modeling. 

D.  MATHEMATICAL  MODELS 

The  general  formulations  of  the  mathematical  models  below  are  developed  in 
previous  works  outlined  in  the  respective  subsection’s  reference  remarks.  The 
presentation  of  the  equations  represents  general  summaries  of  their  work. 

1.  Bilinear  Element 

The  bilinear  element’s  hysteretic  loop  is  shown  above,  in  Figure  5,  and  is  an 
approximation  to  the  dynamic  hysteretic  loop  as  shown  in  Figure  4  above.  The  hysteretic 
loop  is  a  function  of  displacement  and  velocity  across  the  isolator,  primary  spring 
stiffness,  and  the  secondary  spring  stiffness.  Program  code  for  the  hysteretic  loop  is 
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formulated  using  basic  logical  algorithms,  and  is  executed  at  the  given  time  step  interval 
over  the  time  span. 

In  developing  the  algorithm  for  the  hysteresis  loop  it  is  evident  that  the  loop  is  a 
combination  of  4  intersecting  linear  lines.  Additionally,  the  restorative  force  is  a  function 
of  the  spring  constant,  displacement,  velocity,  previous  deformation  state  position,  and 
the  intercept.  Ultimately,  to  identify  which  part  of  the  hysteresis  loop  that  the  current 
deformation  data  is  on,  a  set  of  memorized  variables  is  required.  The  memorized  set  of 
variables  stores  which  leg  the  last  deformation  state  was  at  and  the  force  intercepts  of  the 
4  lines. 

As  an  example,  an  isolator  in  equilibrium  is  subjected  to  a  sinusoidal  base 
excitation.  Assuming  that  the  excitation  begins  in  the  positive  direction  the  initial 
velocity  and  displacement  are  positive.  Given  the  primary  stiffness  as  Ki,  the  force  is 
calculated  from  Equation  (9). 

F=K]Ax  +  Bl 
(9) 

Since  the  equilibrium  position  was  at  zero  velocity  and  zero  displacement  the  value  of  Bi 
is  zero.  If  the  calculated  force  from  Equation  (9)  is  less  than  a  specified  positive 
transition  force,  Fmax,  the  intercept,  B2,  for  the  subsequent  leg  is  simultaneously 
calculated.  This  is  conducted  by  numerically  solving  for  the  linear  equation  with  a  slope 
of  Ko  that  intersects  the  calculated  F  and  Ax  position.  When  the  force,  calculated  in 
Equation  (9),  is  greater  than  the  transition  force,  Fmax,  Equation  (10)  is  used. 

F=K2Ax+B2  (10) 
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Continuing  along  this  leg  the  slope  of  the  next  leg,  B3,  is  simultaneously 
calculated  until  the  velocity  becomes  negative  and  the  isolator  begins  to  unload.  At 
which  point  Equation  (1 1)  is  used. 

F=-KlAx+Bi  (11) 

As  the  deformation  progresses  along  this  leg  the  next  leg’s  intercept,  B4  is 
calculated  as  before.  Upon  surpassing  the  negative  transition  force,  -Fmax,  Equation  (12) 
is  used. 

F  =  -K2Ax+B4  (12) 

The  loop  continues  until  the  end  of  the  time  history  is  reached.  Additionally,  at 
each  time  step,  the  intercepts  and  a  location  variable  are  updated  and  stored. 
Furthermore,  reloading  and  unloading  along  the  same  leg  is  also  possible  given  the  data 
stored  at  each  step. 

The  two  distinct  types  of  bilinear  elements  used  in  this  paper  are  referred  to  as 
“ideal”  and  “real.”  The  real  isolator  models  the  basic  stress-strain  curve,  where  the 
isolator  is  allowed  to  strain  harden.  Therefore,  upon  successive  reload  in  either  the 
tensile  or  compressive  direction  the  yield  force  is  increased.  The  ideal  isolator  is  not 
allowed  to  strain  harden  and  hence  has  constant  positive  and  negative  transition  forces. 

The  algorithm  of  the  example  above  is  presented  in  Figure  6.  Incorporating  this 
algorithm  into  the  recursive  block-by-block  formulation  is  easily  conducted  because  the 
stored  data  is  passed  from  block  to  block.  However,  since  the  process  within  a  specific 
block  is  iterative  the  final  block’s  data  is  not  updated  until  the  respective  block  has 
converged. 
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Figure  6.  Bilinear  Algorithm 
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2.  Wen  Element 


As  stated  earlier,  the  Wen  element  is  a  mathematical  equation  that  accurately 
models  the  hysteretic  behavior  of  nonlinear  isolators.  The  Wen  nonlinear  equation  takes 
the  bilinear  hysteretic  loop  from  Figure  4  and  generates  a  more  realistic  hysteresis.  The 
Wen  element  refines  the  bilinear  element’s  hysteretic  loop  through  variation  of  the  input 
parameters,  as  seen  in  Figure  7. 


Figure  7.  Parameters  of  the  Wen's  Model  From  Ref.  [14:p.  355] 

The  non-dimensional  coefficients  A,  p,  £,  and  T|  control  the  shape  of  the 
hysteresis  loop.  Parameter  A  influences  the  initial  loading  slope  and  subsequent 
unloading  and  reloading  slopes.  The  ratio  of  £/p  controls  the  unloading  and  reloading 
curvatures.  The  effect  of  C/P>1  is  equivalent  to  increasing  the  (A)  parameter  and  causing 
the  loop  to  curve  outward  as  seen  in  Figure  7  where  £>1.  The  parameter  T|  controls  the 
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curves  that  join  the  load  and  unload  curves  and  is  best  described  as  the  elastic-plastic 
behavior  of  the  material. 

The  general  formula  for  calculating  the  restoring  force  is: 

F 

F(t)  =  a—  x(t)  +  (1  +a)Fyz(t)  (13) 

F 

Where  a~r  S(t)  is  the  equivalent  linear  portion  and  (1  +a)FZ(t)  is  the 

o  > 

r 

nonlinear  portion.  The  symbol  a,  is  the  respective  yield  to  post  ratio  parameter,  and  z(t) 
is  the  dimensionless  hysteretic  displacement  that  is  defined  by  the  nonlinear  equation: 


Where,  x  (t)  represents  the  strain.  Rearranging  Equation  (14)  and  substituting  the 
variable  notation  of  £  with  yin  Equation  (14)  the  equivalent  Equations  (15)  and  (16)  are 
adopted  from  [Ref.  13:p.  250]. 

z(0  =  -}|x(t)|z(0"~llz(0l-  fix(tjz(tf  +AS(t )  for  n  =  even  (15) 

z(  t)  =  -  )j x(t)  z(t) n  -  fix  (OkOl"  +  A  x  (0  for  n  =  odd  ( 1 6) 

A  more  refined  formula  for  the  dimensionless  hysteretic  displacement  z(t),  from 
[Ref.  15:p.  79],  is  used  in  this  study. 

z=-7jxjddn-1  -fix\z\n+Ax  (17) 

The  combination  of  the  nonlinear  parameters  from  Equation  (17)  yields  the 
maximum  value  of  the  restoring  force. 
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To  reinforce  the  effects  of  each  of  the  respective  Wen  nonlinear  equation 
parameters  graphical  representations  are  presented.  An  arbitrary  hysteretic  isolator  was 
chosen  and  a  sinusoidal  displacement  was  imposed.  Figure  8  illustrates  the  hysteresis 
loop  for  the  isolator  and  sinusoidal  excitation. 


Figure  8.  Wen  Element  Hysteresis  Loop  for  Fy=l,  A=2,  t|=2,  (3=1,  y=l,  6y=l,  a=.01  due 

to  Unit  Sinusoidal  Excitation  at  .5  Hz 
By  varying  the  constant  parameters  in  Equations  (14)  and  (17),  the  effects, 
delineated  by  the  dashed  lines,  are  presented  in  Figure  9. 
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Figure  9.  Parameter  Change  Effects  On  Wen  Hysteresis,  a)  Increase  Fv  b)  Increase 
A,  c)  Increase  (f  d)  Increase  rj,  e)  Increase  6y  f)  Increase  a ,  g)  Increase  y 
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Maxwell  Element 


The  Maxwell  element  is  another  type  of  nonlinear  isolator  that  is  constructed 
from  a  linear  spring  and  viscous  damper  in  series  vice  in  parallel  as  in  a  common  shock 
absorber.  The  formulation  was  developed  in  [Ref.  14]  and  is  included  in  the  code  to 
provide  an  additional  nonlinear  element.  Although  viscous  dampers  and  shock  absorbers 
are  frequency  dependant,  the  overall  nonlinearity  makes  them  sufficiently  more  efficient 
in  energy  dissipation  then  a  linear  spring.  By  placing  the  spring  and  viscous  damper  in 
series  the  isolator  becomes  frequency  independent. 

The  nonlinear  restorative  force  equation  is: 

F(t)=--F(t)  +  kS(t)  (19) 

c 

E.  LOAD  DEPENDANT  ISOLATORS 

Devices  that  are  used  in  structural  applications  not  only  require  horizontal 
stiffness  but  also  vertical  stiffness  as  well.  The  initial  loading  as  well  as  uplift  loads  on 
the  isolator  plays  an  important  role  in  establishing  the  horizontal  or  shear  stiffness  of  the 
bearing.  The  broad  effects  of  vertical  displacements  on  rubber  and  elastomeric  isolators 
are  well  documented  in  [Ref.  9:p.  170],  [Ref.  16],  and  [Ref.  17]  and  are  broadly 
summarized  in  this  section. 

For  an  isotropic  material  the  loading  versus  the  stiffness  function  is  similar  to  the 
curve  outlined  in  Figure  1 1 . 
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Figure  1 1.  Isolator’s  Properties  Versus  Vertical  Displacements 
For  small  vertical  displacements  the  change  in  the  isolator  loads  will  be  small 
enough  that  the  parameters  will  be  dictated  by  it’s  elastic  range.  Since  the  vertical  spring 
component  of  the  isolator  will  generally  obeys  Hooks  law. 


a=Es  (20) 

Where  c?  represents  the  stress,  s  the  strain,  and  E  the  elastic  modulus.  The  properties  of 
the  isolator  are  then  determined  by  the  vertical  displacements. 

Generally,  when  the  isolator  has  been  loaded  it  has  stored  elastic  energy.  Rubber 
is  particularly  strain  sensitive  which  causes  the  shear  modulus,  G,  to  change  according  to 
the  strains  imparted  due  to  the  connections  to  the  building  and  ground  bearing, 
represented  by  Figure  12. 
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Figure  12.  Shear  Stresses  on  a  Load  Isolator  From  [Ref.  9:p.  142] 

The  comparison  of  the  total  elastic  stored  energy  caused  by  the  shear  forces  on 
the  top  and  bottom  surfaces  to  the  energy  obtained  from  the  average  shear  strain  within 
the  isolator  will  influence  the  overall  effective  shear  modulus. 

rl,=-^Jt\\yl+rl)dydA  (21) 

Where,  A  is  the  cross  sectional  area  of  the  isolator,  and  t  the  thickness.  Depending  on  the 
type  of  isolator,  the  values  of  the  average  shear  strain  can  be  calculated  and  the  resultant 
G  could  then  be  iterated  since  G  is  also  a  function  of  sc  (compression  strain). 

Expanding  the  system  to  more  complex  systems  such  as  anisotropic  materials 
like  composites  requires  a  more  detailed  knowledge  of  the  property  versus  loading 
curves.  Even  a  simple  model  can  lead  to  complex  numerical  formulations  if  the 
transition  was  made  into  the  plastic  region.  This  could  easily  lead  to  a  hysteretic  behavior 
in  the  vertical  component  due  to  the  continuous  loading  and  reloading  in  the  plastic 
range. 
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III.  EQUATIONS  OF  MOTION  AND  NONLINEAR  ISOLATORS 


Generating  the  equations  of  motion  of  an  arbitrary  n-degree  of  freedom  structure 
allows  the  structure  to  be  described  both  statically  and  dynamically  in  a  finite 
mathematical  manner.  By  modeling  the  structure  and  directly  controlling  the  number  of 
degrees  of  freedom  (DOF),  the  number  of  differential  equations  becomes  finite  and 
manageable.  The  left  hand  side  of  the  equations  of  motion  becomes  the  representation  of 
the  elastic-linear  structure  while  the  right  hand  side  captures  all  the  forces  generated 
from  the  linear  and  nonlinear  attributes  of  the  isolators. 

A.  GOVERNING  EQUATIONS 

Developing  a  FE  model  of  a  structure  is  covered  in  various  texts  [Ref.  18]  and 
the  basic  structural  dynamic  formulation  is  covered  in  [Ref.  19].  The  work  in  this  section 
is  a  basic  summary  of  work  developed  in  these  references. 

Taking  a  continuous  model,  which  has  an  infinite  number  of  DOF,  and 
converting  it  to  a  workable  yet  realistic  mathematical  model  is  the  goal  of  FE  techniques. 
In  order  to  determine  the  discretization  of  the  structure  a  logical  determination  has  to  be 
made  in  respect  to  retaining  a  sufficient  number  of  degree  of  freedoms.  The  structure  can 
then  be  describe  by  a  continuous  mathematical  model  and  is  formulated  as  a  partial 
differential  equation  (PDE).  Avoiding  the  computational  effort  associated  with  a  PDE, 
the  system  is  transformed  into  a  set  of  2nd  order  ordinary  differential  equations  (ODE). 
The  transformation  is  conducted  by  using  elements  of  lumped-parameter  models.  The 
formulation  is  covered  in  many  vibration  and  dynamic  analysis  textbooks  as  in  [Ref.  6] 
and  [Ref.  19]  and  is  summarized  below.  Additionally,  it  is  computationally  useful  to 
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transform  the  coupled  ODE  into  modal  space  to  achieve  uncoupled  ODEs.  Again,  these 
methods  are  well  documented  and  are  used  in  summary  in  this  formulation. 

For  an  arbitrary  base  isolated  structure  the  free-free  structure  is  connected  to  the  base 
isolators  as  shown  in  Figure  12. 


n‘ floor 


'  floor 


3re  floor 


2nd  floor 


i  1  floor 


i -^^zzzzz^zzzzz^zbzzz^sbzz  T^^zzzzsLbzzzz^zzzzzJ^za 
Figure  12.  N-Story  Base  Isolated  Structure 

The  structure’s  nodal  degrees  of  freedom  that  connect  to  the  isolator  are 


contained  in  the  vector  (Cset).  (Cset)  can  contain  all  six  degrees  of  freedom  per  node  but 
generally  the  lateral  and  vertical  degree  of  freedoms  are  chosen  for  base  isolation. 


The  governing  equations  of  motion,  assuming  proportional  structural  damping, 
[Ref.4:p.  176],  for  the  coupled  free-free  structure  and  the  isolators,  is  presented  in 
Equation  (22): 

[ML.pl  +ICL.P1  +[*L{rL=feL  (22> 

L  J  nxl  l  Jnxl 

Where,  (n)  represents  the  total  number  of  degrees  of  freedom  in  the  system  and 
the  forcing  function,  F^f°et(t),  contains  the  external  isolator  forces  of  length  (Cset)  and 
zeros  for  non -coupled  degrees  of  freedom.  The  general  solution  to  the  second  order 
differential  equation  is  of  the  form: 

x(t)  =  x(t\+x(t)p  (23) 

The  homogenous  solution  for  particular  degree  of  freedom  becomes: 

x  =  e~^"'  (Asin-Jl-  ^o)nt  +  B  cos-^1  ~Ccont)  (24) 

The  derivation  of  the  particular  solution  depends  upon  the  arbitrary  forcing 
function.  For  simple  forces  such  as  linear  springs  the  solution  becomes  trivial  through 
the  manipulation  of  Equation  (22).  However,  for  nonlinear  elements  the  particular 
solution  can  be  significantly  complex. 

A  method  to  directly  solve  Equation  (22)  is  to  transform  the  second  order 
differential  equation  into  two  coupled  single  ODEs.  Pre-multiplying  Equation  (22)  by 
[M]'1  and  moving  the  stiffness  and  damping  components  to  the  right  yields: 

jxj  =  4M]-'[c]|xJ-[Mr1Mx}+[MrIte,}  (25) 

The  two-coupled  single  ODEs  that  represent  Equation  (25)  are: 
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X 


(26) 


{ j]  =  -WV '  [c]|i|  -  [MX'  [k:Kx}  +  [MX'  te  ] 

The  solution  to  Equation  (26)  can  be  easily  completed  using  any  validated  ODE 
solver,  however,  it  can  become  significantly  problematic  if  the  forcing  matrix 
IXT  {C* }  is  nonlinear. 

B.  VOLTERRA  INTEGRAL  FORMULATION 

The  solution  of  Equation  (22)  can  also  be  obtained  using  the  Volterra  integral 
method.  The  method  of  solving  for  the  particular  solution  of  Equation  (22)  is  conducted 
using  the  convolution  integral.  The  derivation  of  the  convolution  integral  is  developed  in 
various  vibration  textbooks  such  as  [Ref.  6],  and  is  formulated  as: 

x  (t)  =  £  /(  r)h(t  -  z)d t  (27) 

Where,  h(t-x),  is  the  impulse  response  function  defined  as: 

KO  =  ^~=  - -sin  mn  yjl-£2t  (28) 

F 

Where  con  represents  the  natural  frequency,  and  C  the  damping  coefficient.  In 
order  to  sum  each  discretized  time  band  to  yield  the  entire  displacement  time  history  the 
system  must  be  physically  realizable: 

x(t)  =  Xf(r)Kt  -t)At=  lf(j)h(t  -  r)dr  (29) 

Substituting  the  particular  solution  into  Equation  (23),  the  complete  solution  to 
the  general  equation  of  motion  reduces  to: 
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hn(t-T) 

hx,(t-  r) 

-  K(t-r) 

hi2(t-r) 

h22(t-r) 

A„(*-*) 

KSt-t) 

-  hjt-r\ 

F/W 

F/(t) 


\dz 


(30) 


The  resultant  equation  is  the  non-homogenous  Volterra  integral  equation  as  the 
forcing  vector  F/’(t)  is  a  (nonlinear)  function  of  the  response  and  its  derivative.  The 
forcing  vector  F/°(r)  represents  the  isolator  force  due  to  ground  excitation  on  the  DOF 
denoted  in  the  subscript.  In  order  to  reduce  the  computational  effort  for  Equation  (30), 
Equation  (22)  is  transformed  into  modal  space.  The  transformation  equation  is  defined 


as: 


;'= i 


(31) 


Where,  n  is  the  number  of  modes,  i  is  the  i*  DOF,  and  j  is  the  j111  mode. 

Pre-multiply  Equation  (24)  with  the  transpose  of  the  mass  normalized  mode 

T 

shape  matrix,  [<3>]  ,  the  modal  space  equivalent  is: 


(32) 


The  modal  solution  of  Equation  (32)  yields  the  equivalent  Volterra  integral  form 
of  Equation  (33): 
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Since,  the  system  of  equations  was  transformed  into  modal  space  the  impulse 
response  matrix  is  uncoupled  and  is  denoted  by  the  (~)  symbol.  Transforming  Equation 
(33)  back  into  normal  space  by  substituting  in  Equation  (31)  yields: 


Equation  (34)  contains  (n  +  m)  unknowns  and  (n)  coupled  equations,  where  (n) 
represents  the  total  number  of  nodes  and  (m)  the  number  of  nodal  values  contained  in 
(Cset).  Clearly,  we  have  more  unknowns  than  equations  to  easily  substitute  values  in  and 
achieve  a  simple  solution.  However,  since  there  are  only  (m)  linear  or  nonlinear  forces,  a 
recursive  iterative  solution  process  can  be  conducted  solely  on  a  reduced  (m)  sized 
system  of  coupled  equations.  The  set  of  reduced  nodes,  (Cset),  contains  the  associated 
connected  base  nodes  denoted  on  the  subscript  of  F^°(t).  Therefore,  the  number  of 
equations  reduces  to  (m)  and  the  number  of  unknowns  reduces  to  (2m).  This  reduction 
can  dramatically  reduce  the  solution  time  in  comparison  to  methods  that  solve  the  entire 
(n)  governing  ordinary  differential  equations. 


Therefore,  Equation  (34)  reduces  to: 


By  employing  modal  participation  methods,  the  computational  effort  in  solving 
Equation  (57)  can  further  be  reduced.  It  is  shown,  [Ref.  4:ch.  6],  that  for  a  system 
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transformed  into  the  frequency  domain  that  for  a  specified  nodal  maximum  displacement 
to  a  particular  (p)  mode  is  as  follows: 


[x  | 

i  /I  max 


(36) 


Where,  the  superscript  identifies  the  mode,  the  subscript  identifies  the  node,  and  (n)  the 
total  number  of  modes.  From  the  resultant  frequency  response  spectrum,  the 
contributions  of  the  modes  closest  to  the  mode  in  question  will  have  the  most  significant 
impact.  Therefore,  the  transformation  matrix  [<X>],  is  further  reduced  from  a  (n  x  m) 


matrix  into  one  that  is  (m  x  r).  Where,  (r)  represents  the  number  of  modes  retained. 

The  time  domain  normal  mode  truncation  method  approaches  the  number  of 
modes  retained  in  a  similar  manner.  For  example,  using  the  i*  DOF’s  modal 
displacement  time  history  solution,  defined  by  Equation  (39),  it  is  obvious  that  the  more 
modes  retained  the  accuracy  of  the  solution  increases. 

xt{t)  =  O’tf,  (0  +  O,  q2  (?)  +  <%,  (?)  +  ...O; qn  (?)  (37) 

However,  truncating  to  a  sufficient  number  of  modes  with  an  acceptable  error,  can 
reduce  the  computational  effort  dramatically. 

Equation  (35)  now  reduces  to  the  reduced  Volterra  integral  form: 

(r) 

(T)  </r(38) 

(r) 

In  order  to  reduce  Equation  (38)  to  represent  a  single  i*11  DOF  time  history,  the 
equation  is  simplified  by  multiplying  the  matrices  [d>][h][  0]T.  Multiplying  out  [3>][h] 


yields: 
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Continuing,  the  time  history  nodal  response  for  the  ith  DOF  yields: 


xi(t)=xih(t)  + 


f* t 

hjj«- ^so(t)dr+j‘±<fy2'hjj(,- T)F‘so(T)dr+... 
J  J= I  j= 1 


C‘  r 


(40) 

J  0 

Where,  the  notational  superscript  indicates  the  respective  retained  modes  and  the 
subscript  identifies  with  the  DOF  values,  (i),  in  (Cset).  Therefore,  the  summation 
notation  for  (j)  signifies  that  the  associated  DOFs,  contained  in  (Cset),  is  summed  over 
the  total  number  of  modes  retained  from  truncation. 


C.  RECURSIVE  ITERATION  FORMULATION 


By  employing  modal  truncation  and  retaining  the  (Cset)  sized  set  of  equations, 
the  overall  computational  effort  has  been  dramatically  reduced.  The  solution  of  the  (un¬ 
coupled  equations  is  conducted  in  a  recursive  iterative  manner.  From  Equation  (40),  the 
transient  response  for  a  particular  nodal  degree  of  freedom  is  dominated  by  the 

convolution  integrals.  The  isolator  forces,  within  the  respective  integrals  are 

unknown.  Therefore,  in  order  to  solve  for  the  x,(t)  nodal  degree  of  freedom  time  history 
an  initial  guess  is  conducted.  For  a  simple  linear  isolator  the  linear  approximation  can  be 
represented  as: 

Fr=-Kiso{Xi-xs)  (41) 
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Where,  Kiso  represents  the  linear  stiffness,  the  subscript  (i)  the  nodal  DOF,  and  xg 
the  ground  motion.  Substituting  all  the  isolator  forces  into  Equation  (40)  yields: 


*•(0  = 


r  .  •  ~ 

h  jjit-TK-K^ix ,(T) -Xg(T))dT+... 

J  0  j=] 


r  .  .  ~ 

X  ®]m  hjj(f-  r)(-Kiso  Mr)  -  xs(T)))dT  (42) 

J  o  M 

The  solution  for  any  one  nodal  DOF  then  becomes  a  function  of  all  the  (Cset) 
unknown  DOF  displacements.  To  achieve  a  solution  a  recursive  interative  process  is 
employed  for  each  nodal  DOF.  Therefore,  Equation  (42)  is  rewritten  as: 


*/ 


'(0  = 


r  r 


o;^1 


h jjit-  r)(-KrJx“' (T) - xt(T))dT+... 


W  >  -  *)(-**(*“(*)  -  (.mdr  (43) 

J  0  J=l 

The  recursive  iteration  process  uses  the  x^s'(z)  values  to  solve  for  the  xn™ (?) 
values.  The  iterative  process  begins  by  assigning  initial  arbitrary  values  to  x^\t) 
through  x‘“s'(z).  These  values  are  then  substituted  into  Equation  (43)  and  the  resultant 
x™w(z)  values  are  computed.  A  tolerance  check  is  conducted  between  x!“st(z)  and 
x”ew(' z),  if  the  criteria  is  not  met  the  M(z)  values  are  updated  to  the  x™w(x)  values. 
The  process  continues  until  tolerance  is  met.  An  algorithm  of  the  iterative  recursive 
method  is  presented  in  Figure  13. 
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End 

Figure  13.  Basic  Recursive  Iteration  Algorithm 

D.  LINEAR  AND  NONLINEAR  RECURSIVE  BLOCK  BY 
BLOCK  SYNTHESIS 

The  algorithm  used  to  solve  the  governing  Volterra  integral  can  still  be 
computationally  exhaustive  for  large,  finely  discretized  time  spans  and  or  large  (Cset) 
DOF  systems.  For  physically  realizable  systems  the  computational  effort  can  be  further 
decreased  by  recursive  block-by-block  iteration.  The  convolution  integral  can  be  divided 
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into  blocks  and  summed  up  to  arrive  at  the  total  solution.  The  blocks  take  advantage  of 
the  linear  properties  of  the  convolution  integral. 

By  the  formulation  of  the  convolution  integral,  Equation  (27),  the  transient 
response  spanned  over  some  period  of  time  t,  is  the  summation  of  the  time  histories  of 
the  arbitrary  impulse  forces  acting  at  specified  times,  x. 

~  *)dr  =  f(0)h{t -  r)dr  +  |  /(/,  )h(t  -  z)dz  + ... |  / (tn)h(t  -  r)dz  (44) 

Therefore,  the  incorporation  of  this  linearity  of  the  convolution  integral  into  the 
recursive  block-by-block  method  is  presented.  In  general,  each  block  will  synthesize 
both  a  specified  transient  time  history  and  an  associated  forcing  function.  The  recursive 
process  requires  the  use  of  all  previous  block’s  solutions  to  achieve  the  total  solution.  To 
represent  the  process  Equation  (44)  can  be  rewritten  to  represent  discrete  block  times  as: 

| \f{r)Kt-r)dr=  ^f(0...tl)h(t-r)dr  +  ^f(tv.i2)h(t -  r)dr+ ... 

[  f(tn...t)h{t -t)dr  (45) 

Where,  the  forcing  functions  /(0.../,)  represent  the  block’s  synthesized  force 
over  the  specified  block  time  span.  The  recursive  properties  are  not  outwardly  evident  in 
the  formulation  of  Equation  (45)  but  is  accounted  for  computationally  for  each 
successive  block.  To  illustrate  the  process  a  system  is  synthesized  using  two  blocks. 
Equation  (45)  reduces  to: 

{f(TM  ~  T)dr  =  -  r)dr  +  |  fit,  ..I)h(f  -  r)dr  (46) 
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Computationally,  Equation  (44)  is  iterated  and  converged  for  the  first  block.  The 
next  block  then  proceeds  to  solve  Equation  (44)  for  the  second  time  block.  However,  the 
first  block’s  isolator  forcing  time  history  solution  has  to  be  extended  over  the  second 
block: 

Here,  h(tv..t-  r)  represents  the  impulse  response  function  that  begins  at  time  fi 

vice  zero.  The  time  history,  x(t),  for  the  extended  time  block  J  z)dr  is 

the  of  the  next  block.  Although  it  can  be  generally  said  that  any  arbitrary  x^'(r) 

can  be  chosen  for  each  successive  blocks,  a  closer  approximation  to  the  solution  can 
reduce  the  iterative  process.  More  importantly,  the  computed  extension  is  added  to  the 
final  recursive  iterated  solution  of  the  second  block.  The  recursive  block-by-block 
method  algorithm  is  presented  in  Figure  14. 


of(r)h(t-T)dr=  ^  f(0...tA)h(t -  r)dz+  )h(tx...t - r)d.T  (47) 
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End 


Figure  14.  Recursive  Block-by-Block  Algorithm 
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For  linear  isolators  such  as  the  one  presented  in  Equation  (41)  block-by-block 
iterations  and  block  summations  are  straightforward.  Each  block  follows  the  algorithm 
presented  in  Figure  16  and  Figure  17.  However,  for  nonlinear  forces  the  recursive  block- 
by-block  synthesis  relies  more  heavily  upon  the  force  solutions  of  the  previous  blocks. 

A  possible  model  of  a  nonlinear  isolator  is  presented  as: 

F/'0  =f((xg-x]),(xg-x’),(xg-x.),Ki(t),....)  (48) 

Where  the  superscript  *  symbol  indicates  the  iterated  value  achieved  from  the 
algorithm  presented  in  Figure  16.  For  the  nonlinear  isolators  presented  a  Runge-Kutta  4, 
(RK4),  ODE  solver  is  used. 

Each  successive  block  requires  that  the  solution  to  the  nonlinear  equations  have 
the  correct  initial  conditions.  Therefore,  in  order  to  achieve  continuity  of  the  isolator’s 
hysteretic  response  between  blocks,  the  initial  condition  is  interpolated  through  a  direct 
quadratic  approximation  through  the  last  three  force  data  points  of  the  previous  block,  as 
shown  in  Figure  15. 
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Figure  15.  Recursive  Interpolation 


E.  UPLIFT  IN  RECURSIVE  BLOCK-BY-BLOCK  SYNTHESIS 


Beginning  with  Equation  (42)  the  calculation  of  uplift  requires  the  retention  in 
the  analysis  of  the  vertical  degrees  of  freedom.  Therefore,  if  the  original  retained  set  of 
nodal  degrees  of  freedom  were  limited  to  lateral  displacements  in  (Cset),  the  new  (Cset) 
would  include  the  vertical  DOFs  associated  with  the  Cset  nodes.  This  would  increase 
(Cset)  by  a  factor  of  two.  Including  bi-directional  excitation  with  uplift  would  then  result 
in  tripling  the  size  of  (Cset)  [Ref.  20]. 

For  different  isolators  the  uplift  properties  and  their  effects  on  lateral  stiffness  is 
well  documented  in  [Ref.  21],  [Ref.  16],  and  [Ref.  17].  For  simplicity  the  uplift  portion 
of  the  synthesis  was  modeled  as  a  linear  spring: 

r —KiS0(Xrcal)  (49) 

Programming  interpolating  tables  can  easily  be  introduced  in  the  developed 
program  if  known.  Therefore,  Equation  (40)  expands  to  include  uplift  synthesis: 

x.  (/)  =  Equation^  42)  +  1 2>/< + !  ~hu  <'  -  *Kpl + - 

7=1 
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Additionally,  since  many  of  the  parameters  are  uplift  dependant  the  following 


isolator  models  are  presented: 


/r=-c  (*,-*,) 


(51) 


F’s°  =  /((x  -x*),(xg-  x'),(xg-  x'),K'(t),  j]* 


(52) 


Where,  the  superscript  symbol,  *,  represents  a  now  synthesized  value  as  depicted 
in  the  linear  approximation  of  slope  (M)  below: 

Paramater-  (t)  =  Parameter™*™1  +  M(xJer,ical( 0)  -  xjer,ica! (t))  (53) 

F.  SINGLE  ISOLATION  TRAPEZIODAL  EXTRACTION 


Given  an  arbitrary  excitation  displacement  time  history  the  resulting  forcing 
function  can  be  extracted  using  a  trapezoidal  integration  scheme.  The  general  formula 


taken  from  any  numerous  numerical  texts  [Ref.  22]  is  presented  as: 

x(t)=— A t{fx  +  2f2+2f3+  2 /4...  +  /„)  (54) 

Applying  the  formulation  to  the  convolution  integral  yields: 

x(0)  =  ^AP(0)/(0))  (55) 

x(A t)  =  |  At(/<0)/(0)  +  KAt)f(At))  (56) 

x(nAt)  =  -^dt(h(0)f(0)  +  2h(At)f(At)  +  ...h(nAt)f(nAt)')  (57) 


Rearranging  Equations  (55)  through  (57)  into  matrix  form  yields: 


(58) 


'  *(0)  ' 

'  0 

0 

•  •  • 

0 

O' 

'  /(O) ' 

x(At) 

■  =  -Al 
o 

h(At) 

0 

•  •  • 

0 

0 

/(At) 

x(2At) 

h{  2  At) 

h(At) 

•  •  • 

0 

0 

f(2At)  > 

• 

jL 

• 

• 

•  •  • 

•  ••  0 

• 

x(nAt ) 

h(nAt ) 

h((n-  l)At) 

•  •  • 

h{At)  0_ 

f(nAt) 

The  impulse  response  function  matrix  is  a  lower  triangular  matrix  that  has  zeros 
dominating  the  main  diagonal.  This  in  effect  negates  the  weighting  factor  of  the  1/2.  The 
integration  scheme  then  turns  into  a  rectangular  integration  method  with  an  0(h)  global 
error.  To  achieve  the  same  0(h  )  approximation  of  the  trapezoidal  scheme  the  weight  is 
placed  on  f(0).  Rearranging  the  matrices  in  order  to  solve  for  f(t)  yields: 


'  x(A 0 ' 

'  h(At) 

0 

•  ••  0 

'  2/(0)  ' 

x(2At) 

h(At) 

h(At) 

•  ••  0 

/(At) 

x(3At) 

>=At 

h(At) 

0 

h(At) 

0 

•  ••  o 

0  0  0  0 

■  /(2At) 

0 

x(nAt) 

h(nAt ) 

h((n  -  l)At) 

•  •  •  h(At) 

/((n-l)At) 

The  solution  for  the  extracted  forcing  function  can  be  accomplished  through 
various  numerical  methods.  The  trapezoidal  integration  scheme  is  a  simple  way  of 
acquiring  the  synthesized  forces  given  a  non-isolated  simple  structure.  Increasing  the 
time  step  would  increase  accuracy,  but  can  be  computationally  exhaustive  for  large  time 
histories. 
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V.  RESULTS 


To  develop  MATLAB  code  for  the  RBBIEF  method  workable  and  especially 
comparative  models  were  used.  A  sinusoidal  wave  and  the  El  Centro  North  South 
Earthquake  of  1940  were  used  for  the  base  excitation. 

Appendix  A  presents  the  validation  of  the  trapezoidal  method  versus  NASTRAN 
and  using  MATLAB’ s  ode45  function.  Using  the  same  model  the  recursive  block-by¬ 
block  method  is  validated  against  ode45  using  a  simple  spring.  The  model  is  then  run  for 
an  ideal  bilinear  element  and  compared  to  the  results  of  the  recursive  block-by-block 
method 

In  order  to  demonstrate  the  advantages  and  disadvantages  of  the  recursive  block- 
by-block  method  3  models  were  synthesized  using  the  1940  El  Centro  earthquake  data. 
The  first  finite  element  model  is  a  4-story,  single  square  bay,  that  is  25  feet  wide  and 
having  an  inter-story  height  of  17.5  feet.  The  structure  is  modeled  using  a  total  of  20 
nodes,  120  DOF,  four  DOF  in  excitation,  and  four  DOF  in  uplift.  The  second  finite 
element  model  is  a  4-story,  nine  square  bays  that  is  have  a  total  width  of  75  and  having 
an  inter-story  height  of  17.5  feet.  The  building  is  modeled  using  a  total  of  80  nodes,  480 
DOF,  with  16  DOF  in  excitation,  and  16  DOF  in  uplift.  The  last  model  is  the  30-story 
structure  that  is  225  feet  by  225  feet  at  the  base  and  an  inter-story  height  of  17.5  feet  for 
a  total  of  525  feet  in  height.  The  30-story  structure  is  modeled  using  2,51 1  nodes,  15,066 
degrees  of  freedom,  81  DOFs  in  excitation,  and  81  DOFs  in  uplift.  A  summary  is 
presented  in  Table  1. 
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1  -Bay  4-Story 

9-Bay  4-Story 

30-Story 

Nodes 

20 

80 

2511 

Deqrees  of  Freedom 

120 

480 

15066 

DOF  in  exictation 

4 

16 

81 

DOF  in  Uplift 

4 

16 

81 

Base  Dimension 

25'  x  25’ 

75'  x  75' 

225*  x  225* 

Total  Height 

70' 

70’ 

525’ 

Table  1.  Parameters  of  Modeled  Buildings 


D.  PROGRAM  VALIDATION 

1.  De-coupled  Versus  Coupled  System 

A  similar  recursive  process  is  validated  in  [Ref.  24]  and  is  presented  in  Appendix  B  to 
validate  the  trapezoidal  method.  Using  the  trapezoidal  method,  a  9-DOF  spring  mass 
system,  as  shown  in  Figure  16,  is  used  to  present  the  validity  of  uncoupling  a  system. 
The  equivalent  de-coupled  system  is  illustrated  in  Figure  17.  Using  a  single  block  in  the 
recursive  block-by-block  method  with  a  linear  spring  of  10  lbf/in  as  the  isolator  the 
comparison  is  conducted.  Both  programs  were  run  using  a  time  step  of  .01  seconds  for  a 
ten  second  duration  with  a  sinusoidal  ground  excitation  of  unit  amplitude,  oscillating  at  a 
frequency  of  .5  Hz. 
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Figure  16.  9-DOF  Spring  Mass  System 


Figure  17.  Uncoupled  Large  Mass  9-DOF  Spring  Mass  System 


Synthesizing  a  solution  using  single  block  in  the  recursive  block-by-block 


program,  Figure  18,  represents  the  base  displacement  time  history. 


Nodal  response  of  first  node 


Time  (s) 

Figure  18.  Base  Displacement  Time  History  of  9-DOF  System  Using  Recursive 

Block-by-Block  Synthesis 

The  trapezoidal  results  plotted  along  the  curve  presented  in  Figure  18,  so  absolute 
differences  were  calculated  between  the  two  data  sets  and  revealed  that  the  maximum 
difference  is  on  the  order  of  10'8. 
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2.  Block-by-Block  Synthesis  Versus  MATLAB's  ODE45 


The  recursive  block-by-block  validation  process  will  use  the  simple  9-DOF 
system  is  modeled  in  Figure  16.  The  process  will  begin  with  comparing  the  two  methods 
using  a  simple  linear  spring.  MATLAB’s  ode45  function  will  use  the  state  space 
representation  of  the  second  order  differential  equation  presented  by  Equation  (14).  A 
comparison  with  an  ideal  bilinear  isolator  is  then  made.  Because  ode45  is  an  adaptive 
algorithm,  memory-type  isolators  with  distinct  discontinuities  at  transition  points  along 
the  loop  can  make  the  adaptive  solver  computationally  exhaustive.  For  example,  the 
bilinear  element  has  sharp  transition  points  around  which  the  ode45  function  will  tend  to 
oscillate.  To  avoid  these  oscillations,  additional  code  was  required  to  force  the  solution 
to  stay  within  the  hysteretic  loop. 

Using  the  same  model  in  Figure  16,  an  arbitrary  (nt)  number  of  blocks  were 
chosen  producing  the  same  plot  as  Figure  19.  For  clarification  the  number  of  block  (nt) 
was  set  to  three,  the  linear  spring’s  value  was  10  lbf/in,  and  the  tolerance  in  the  recursive 
block-by-block  synthesis  was  set  at  10'6  of  displacement  and  10"4  for  velocity. 
Conducting  the  ode45  equivalent  and  plotting  the  two  graphs  together  depicted  no 
appreciable  difference.  Therefore  an  absolute  difference  plot  is  again  conducted  in  Figure 
19. 
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Figure  19.  9-DOF  Spring  Mass  System,  Absolute  Base  Displacement  Differences 
Between  a  Three  Block  Recursive  Block-by-Block  Synthesis  and  ODE45 
Continuing  the  comparisons  an  ideal  bilinear  element  was  chosen  next.  The 
maximum  tensile  and  compressive  forces  were  set  to  5  Ibf  and  -5  lbf  respectively.  The 
initial  spring  constant,  Ki,  was  set  to  lOlbf  with  a  yield  to  post  ratio  of  .04.  Again,  the 
ground  excitation  was  a  unit  sinusoidal  waveform  with  a  frequency  of  .5  Hz.  Both 
programs  were  conducted  with  a  time  step  of  .01  seconds  and  a  duration  of  10  seconds. 
Figure  20,  presents  the  two  time  displacement  histories  of  the  bottom  node.  The  dashed 
line  represents  the  ode45  solution.  Figure  21,  represents  the  two  hysteretic  loops.  The 
circles  represent  the  recursive  ode45  solution. 
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Oao  I  ‘lot 


Figure  20.  Three  Block  Synthesis  Versus  ODE45  Base  Displacements 
Comparisons  of  a  Linear  Spring  Isolated  9-DOF  Spring  Mass  System,  from  a  Unit 
Sinusoidal  Ground  Excitation  at  .5  Hz 


Figure  21 .  Three  Block  Synthesis  Versus  ODEA 5  Hysteresis  Comparisons  of  a 
Linear  Spring  Isolated  9-DOF  Spring  Mass  System,  From  A  Unit  Sinusoidal  Ground 

Excitation  At  .5  Hz 
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From  Figure  21  it  is  shown  how  the  ode45  function  can  oscillate  about  the 
discontinuous  transition  comers.  Ultimately,  Figure  20  indicates  a  very  good 
approximation  of  the  synthesis  in  comparison  to  the  ode45  solution. 

3.  Model  Isolator  Plots 

A  sinusoidal  wave  of  unit  amplitude  and  a  frequency  of  .5  Hz  was  used  on  the 
following  modeled  isolators  below  in  Table  2. 


Properties 

Real  Bilinear  Element 

K  =  500  Ibf/in 

Yield  to  post  ratio  =  .04 

Max  Tensile  Force  =  50  Ibf 

Max  Compressive  Force  =  -50  Ibf 

Ideal  Bilinear  Element 

K  =  500  Ibf/in 

Yield  to  post  ratio  =  .04 

Max  Tensile  Force  =  50  Ibf 

Max  Compressive  Force  =  -50  Ibf 

Wen  Element 

Fy=  100  Ibf/in 

Max  Elastic  Displacement  =1  in 

A=8 

n  =  2 
g  =  .0001 
v=.0001 
alpha=.  04 

Maxwell  Element 

K  =  500  Ibf/in 

C  =  51  Ibf-s/in 

Table  2.  Model  Isolator  Paramters  of  Sinusoidal  Input 


Figures  22  through  25  plot  the  respective  hysteretic  response  loops.  The  duration 
of  excitation  was  ten  seconds  and  at  a  time  step  of  .01  seconds.  The  values  in  Table  2  are 
arbitrary  and  were  chosen  to  proved  illustration  in  the  figures. 
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Force  (lbf)  **  Force  (ibf) 
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B.  SINGLE  BAY  4-STORY  BUILDING 


The  1940  El  Centro  Earthquake  Data,  Figure  26,  displacement  time  history  is 
used  for  the  ground  motion  for  synthesis.  The  earthquake  data  was  sampled  at  .1  seconds 
for  a  duration  of  50  seconds. 


Time  (s') 


Figure  26.  1940  El  Centro  North  South  Ground  Motion  Time  History 

The  uplift  isolators  were  modeled  as  linear  springs.  The  columns  and  beams  are  50  ksi 
steel  structural  members  with  the  following  specifications:  columns,  W36x486;  first  and 
second  floor  beams,  W36xl70;  third  floor  beams  are  W36xl60;  fourth  floor  beams, 
W36xl50;  roof  beams,  W36xl35.  A  simple  frame  illustration  of  the  building  is  shown  in 
Figure  27.  The  isolator’s  properties  are  listed  in  Table  3.  Figure’s  28-39 
illustrates  the  hysteresis,  horizontal  and  vertical  displacement  time  responses  to  the 
earthquake  input.  The  direction  of  excitation  is  along  the  weak  axis  of  bending. 
Additionally,  a  cutoff  frequency  of  10  Hz  was  chosen  leading  to  ten  retained  modes. 
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Figure  27.  Single  Bay,  4-Story  Building  Frame 


j|  Properties 

Real  And  Ideal  Bilinear 

Elements 

K  =  300  Ibf/in 

Yield  to  post  ratio  =  .04 

Max  Tensile  Force  =  50  Ibf 

Max  Compressive  Force  =  -50  Ibf 

K  vertical  =  30  Ibf/in 

Stiffness  Slope  =  30  Ibf/in-in 

Wen  Element 

Fy  =  20  Ibf/in 

Max  Elastic  Displacement  =3  in 

A=3 

n  =  4 
g  =.0001 
v=  000001 

K  vertical  =  30  Ibf/in 
alpha=.04 

A,g,n,v,xy,Fy  slope  = 
.08/in*(A,g,n,v,xy,Fy) 
alpha  slope=  008/in*alpha 

Maxwell  Element 

K  =  300  Ibf/in 

C  =  51  Ibf-s/in 

K  vertical  =  30  Ibf/in 

Stiffness  Slope  =  30  Ibf/in-in 

Table  3.  Single  Bay,  4-Story  Isolator  Properties 


Isolator  Displacement  (in) 


Figure  28.  Single  Bay,  4-Story,  Real  Bilinear  Isolator,  Comer  Node  Hysteresis  for  the 
1940  El  Centro  North  South  Earthquake  Ground  Motion 


Lateral  Comer  Node 


Time  (s) 


Figure  29.  Single  Bay,  4-Story,  Real  Bilinear  Isolator,  Comer  Node  Lateral 
Displacement  History  for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 
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Comer  Node  Uplift 


Time  (s) 


Figure  30.  Single  Bay,  4-Story,  Real  Bilinear  Isolator,  Comer  Node  Uplift  Displacement 
History  for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 
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Figure  3 1 .  Single  Bay,  4-Story,  Ideal  Bilinear  Isolator,  Comer  Node  Hysteresis  for  the 
1940  El  Centro  North  South  Earthquake  Ground  Motion 


Lateral  Hysteresis  of  Comer  Node 
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Lateral  Comer  Node 


Time  (si) 


Figure  32.  Single  Bay,  4-Story,  Ideal  Bilinear  Isolator,  Comer  Node  Lateral 
Displacement  History  for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 


Time  (s) 


Figure  33.  Single  Bay,  4-Story,  Ideal  Bilinear  Isolator,  Comer  Node  Uplift  Displacement 
History  for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 
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Figure  34.  Single  Bay,  4-Story,  Wen  Isolator,  Comer  Node  Hysteresis  for  the  1940  El 
Centro  North  South  Earthquake  Ground  Motion 


Lateral  Corner  Node 


Figure  35.  Single  Bay,  4-Story,  Wen  Isolator,  Comer  Node  Lateral  Displacement  History 
for  the  1 940  El  Centro  North  South  Earthquake  Ground  Motion 
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Figure  36.  Single  Bay,  4-Story,  Wen  Isolator,  Comer  Node  Uplift  Displacement  History 
for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 


Isolator  Displacement  (in) 


Figure  37.  Single  Bay,  4-Story,  Maxwell  Isolator,  Comer  Node  Hysteresis  for  the  1940 
El  Centro  North  South  Earthquake  Ground  Motion 
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Lateral  Comer  Node 


Figure  38.  Single  Bay,  4-Story,  Maxwell  Isolator,  Comer  Node  Lateral  Displacement 
History  for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 


Figure  39.  Single  Bay,  4-Story,  Maxwell  Isolator,  Comer  Node  Uplift  Displacement 
History  for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 
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Table  4  and  Table  5,  present  the  computational  results  of  the  synthesis  process 


using  MATLAB  5.2.1  on  a  733  Mhz  Pentium  III,  Dell  Computer. 


Blocks 

Element 

4 

8 

16 

24 

Spring 

8.3 

6 

6.4 

8.8 

Real  Bilinear  Element 

7.6 

5.4 

6.6 

8.6 

Ideal  Bilinear  Element 

5.5 

4.3 

6 

8.5 

Wen  Element 

5.1 

4.6 

6.1 

9 

Maxwell  Element 

8.2 

6.4 

7.25 

10 

Table  4.  Single  Bay,  4-Story,  1940  El  Centro  Earthquake  Recursive  Block-By-Block 

Synthesis  Time  Results  in  Seconds 


Blocks 

Element 

4 

8 

16 

24 

Spring 

5.50E+08 

9.30E+07 

5.17E+07 

4.22E+07 

Real  Bilinear  Element 

1.64E+08 

7.11E+07 

4.78E+07 

3.90E+07 

Ideal  Bilinear  Element 

1.23E+08 

6.00E+07 

4.38E+07 

3.90E+07 

Wen  Element 

2.10E+08 

8.89E+07 

4.40E+07 

3.77E+07 

Maxwell  Element 

1.50E+08 

7.04E+07 

4.93E+07 

4.14E+07 

Table  5.  Single  Bay,  4-Story,  1940  El  Centro  Earthquake  Block-By-Block  Synthesis 

FLOP  Count  Results 

From  the  preceeding  tables  there  is  an  indication  that  although  the  number  of 
FLOPS  decrease  per  block  the  is  an  optimum  number  of  blocks.  This  is  due  to  the 
number  of  optimum  iterations  per  block.  Upon  converging  to  a  number  of  iterations  per 
block  this  number  remains  relatively  the  same  for  each  block.  Therefore,  increasing  the 
number  of  blocks  with  the  same  number  of  iterations  increases  the  computation  time. 
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C.  9-BAY,  4-STORY  BUILDING 


Increasing  the  model’s  complexity,  a  9-bay  4-story  building,  Figure  40,  is  now 
investigated  using  the  16  hysteretic  isolators  that  include  uplift.  Uplift  is  again  calculated 
as  in  the  single  bay,  4-story  structure. 


Figure  40.  9-Bay,  4-Story  Building  Frame 

The  building  is  design  with  50  ksi  steel  members.  The  columns  are  W36x486, 
first  and  second  floor  beams  are  W36xl70,  third  floor  beams  are  W36xl60,  the  fourth 
floor  beams  are  W36xl50,  and  the  roof  beams  are  W36xl35. 


The  structure  is  synthesized  for  various  isolators  listed  in  Table  6,  with  a  base 
excitation  of  the  1940  El  Centro  earthquake.  The  transient  time  history  is  developed 
from  a  time  step  of .  1  seconds  for  a  duration  of  50  seconds. 


Properties 

Real  And  Ideal  Bilinear 

Elements 

K  =  300  Ibf/in 

Yield  to  post  ratio  =  .04 

Max  Tensile  Force  =  50  Ibf 

Max  Compressive  Force  =  -50  Ibf 

K  vertical  =  30  Ibf/in 

Stiffness  Slope  =  30  Ibf/in-in 

Wen  Element 

Fy  =  20  Ibf/in 

Max  Elastic  Displacement  =3  in 

A=3 

n  =  4 
g  =.00001 
v=000001 

K  vertical  =  30  Ibf/in 
alpha=.04 

A.g^v.xy.Fy  slope  = 
.08/in*(Atg,n,v,xy,Fy) 
alpha  slope=.008/in*alpha 

Maxwell  Element 

K  =  300  Ibf/in 

C  =  51  Ibf-s/in 

K  vertical  =  30  Ibf/in 

Stiffness  Slope  =  30  Ibf/in-in 

Table  6.  9-Bay,  4-Story  Isolator  Properties 

Figures  41-52,  illustrate  the  results  of  displacement  time  histories  and  hysteretic 
loops  of  the  respective  isolators.  Again,  the  cutoff  frequency  was  set  to  10  Hz  retaining 
20  modes.  The  direction  of  excitation  is  along  the  weak  axis  of  bending. 
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Figure  41 .  9-Bay,  4-Story,  Real  Bilinear  Isolator,  Comer  Node  Hysteresis  for  the  1940  El 
Centro  North  South  Earthquake  Ground  Motion 


Lateral  Corner  Node 


Figure  42.  9-Bay,  4-Story,  Real  Bilinear  Isolator,  Comer  Node  Lateral  Displacement 
History  for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 
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Figure  43.  9-Bay,  4-Story,  Real  Bilinear  Isolator,  Comer  Node  Uplift  Displacement 
History  for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 


Figure  44.  9-Bay,  4-Story,  Ideal  Bilinear  Isolator,  Comer  Node  Hysteresis  for  the  1940 
El  Centro  North  South  Earthquake  Ground  Motion 
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Lateral  Corner  Node 


Figure  45.  9-Bay,  4-Story,  Ideal  Bilinear  Isolator,  Comer  Node  Lateral  Displacement 
History  for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 


Figure  46.  9-Bay,  4-Stoiy,  Ideal  Bilinear  Isolator,  Comer  Node  Uplift  Displacement 
History  for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 
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Figure  47.  9-Bay,  4-Story,  Wen  Isolator,  Comer  Node  Hysteresis  for  the  1940  El  Centro 

North  South  Earthquake  Ground  Motion 


Lateral  Corner  Node 


Figure  48.  9-Bay,  4-Story,  Wen  Isolator,  Comer  Node  Lateral  Displacement  History  for 
the  1940  El  Centro  North  South  Earthquake  Ground  Motion 
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Corner  Node  Uplift 


Figure  49.  9-Bay,  4-Story,  Wen  Isolator,  Comer  Node  Uplift  Displacement  History  for 
the  1940  El  Centro  North  South  Earthquake  Ground  Motion 


Figure  50.  9-Bay,  4-Story,  Maxwell  Isolator,  Comer  Node  Hysteresis  for  the  1940  El 
Centro  North  South  Earthquake  Ground  Motion 
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Corner  Node  Uplift 


Figure  5 1 .  9-Bay,  4-Story,  Maxwell  Isolator,  Comer  Node  Lateral  Displacement  History 
for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 


Lateral  Corner  Node 


Figure  52.  9-Bay,  4-Story,  Maxwell  Isolator,  Comer  Node  Uplift  Displacement  History 
for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 
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Again,  Table  7  and  Table  8,  present  the  computational  results  of  the  synthesis 


process  using  MATLAB  5.2  on  a  733  Mhz  Pentium  III,  Dell  Computer. 


Blocks 

Element 

2 

4 

8 

16 

Spring 

124 

59 

50 

71 

Real  Bilinear  Element 

50 

49 

44 

47 

Ideal  Bilinear  Element 

61 

38 

39 

65 

Wen  Element  1 

57 

38 

40 

66 

Maxwell  Element 

73 

44 

45 

72 

Table  7.  9-Bay,  4-Story,  1940  El  Centro  Earthquake  Recursive  Block-By-Block 

Synthesis  Time  Results  in  seconds 


Blocks 

Element 

2 

4 

8 

16 

Spring 

9.27E+09 

2.89E+09 

7.37E+09 

1.54E+09 

Real  Bilinear  Element 

3.17E+09 

2.05E+09 

7.01  E+08 

1.45E+09 

Ideal  Bilinear  Element 

4.06E+09 

1.62E+09 

6.74E+08 

1.45E+09 

Wen  Element 

2.70E+09 

1.24E+09 

6.39E+08 

1 .22E+09 

Maxwell  Element 

5.10E+09 

1.69E+09 

6.93E+08 

1.54E+09 

Table  8.  9-Bay,  4-Story,  1940  El  Centro  Earthquake  Block-By-Block  Synthesis  FLOP 


Count  Results 


As  before,  there  indicates  an  optimum  number  of  blocks.  In  comparison  to  the 
single  bay  4-story  structure  the  optimum  number  of  blocks  were  relatively  the  same 
although  the  number  of  FLOPS  increased. 
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D.  30  STORY  BUILDING 


The  30-story  building,  Figure  53,  is  now  investigated  using  81  various  hysteretic 
isolators  that  include  uplift.  Uplift  calculations  are  carried  from  the  previous  examples. 
The  structure  is  synthesized  for  various  isolators  listed  in  Table  9,  with  a  base  excitation 
of  the  1940  El  Centro  earthquake. 


f  Properties 

Real  And 
Ideal  Bilinear 
Elements 

K  =  300  Ibf/in 

Yield  to  post  ratio  =  .04 

Max  Tensile  Force  =  50  Ibf 

Max  Compressive  Force  =  -50  Ibf 

K  vertical  =  30  Ibf/in 

Stiffness  Slope  =  30  Ibf/in-in 

Wen  Element 

Fy  =  20  Ibf/in 

Max  Elastic  Displacement  =3  in 

A=3 

n  =  4 
g  =.00001 
v=.000001 

K  vertical  =  30  Ibf/in 
alpha=.04 

A,g,n,v,xy,Fy  slope  = 
.08/in*(A,g,n,v,xy,Fy) 
alpha  slope=.008/in*alpha 

Maxwell 

Element 

K  =  300  Ibf/in 

C  =  51  Ibf-s/in 

K  vertical  =  30  Ibf/in 

Stiffness  Slope  =  30  Ibf/in-in 

Table  9.  30-Story  Isolator  Properties 

The  transient  time  history  is  developed  from  a  time  step  of  .1  seconds  for  a 
duration  of  50  seconds.  Again,  excitation  is  in  weak  axis  of  bending  and  Figures  54-65 
represents  the  results  of  the  displacement  time  histories  and  hysteresis. 
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Figure  53.  30-Story  Building  Wire  Frame 


Figure  54.  30-Story,  Real  Bilinear  Isolator,  Comer  Node  Hysteresis  for  the  1940  El 
Centro  North  South  Earthquake  Ground  Motion 


Lateral  Corner  Node 


Figure  55.  30-Story,  Real  Bilinear  Isolator,  Comer  Node  Lateral  Displacement  History 
for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 
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Figure  56.  30- Story,  Real  Bilinear  Isolator,  Comer  Node  Uplift  Displacement  History  for 
the  1940  El  Centro  North  South  Earthquake  Ground  Motion 


Figure  57.  30-Story,  Ideal  Bilinear  Isolator,  Comer  Node  Hysteresis  for  the  1940  El 
Centro  North  South  Earthquake  Ground  Motion 
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Lateral  Corner  Node 


Figure  58.  30-Story,  Ideal  Bilinear  Isolator,  Comer  Node  Lateral  Displacement  History 
for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 


Figure  59.  30-Story,  Ideal  Bilinear  Isolator,  Comer  Node  Uplift  Displacement  History 
for  the  1940  El  Centro  North  South  Earthquake  Ground  Motion 
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Lateral  Hysteresis  of  Corner  Node 


-150  1 - 1 - 1 - 1 - 1 - 1 - 1 

-15  -10  -5  0  5  10  15 

Isolator  Displacement  (in) 

Figure  60,  30-Story,  Wen  Isolator,  Comer  Node  Hysteresis  for  the  1940  El  Centro  North 

South  Earthquake  Ground  Motion 


Lateral  Corner  Node 
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0  10  20  30  40  50  60 

Time  (s) 

Figure  61.  30-Story,  Wen  Isolator,  Comer  Node  Lateral  Displacement  History  for  the 
1940  El  Centro  North  South  Earthquake  Ground  Motion 
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Figure  62.  30-Story,  Wen  Isolator,  Comer  Node  Uplift  Displacement  History  for  the 
1940  El  Centro  North  South  Ground  Motion 
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Figure  63,  30-Story,  Maxwell  Isolator,  Comer  Node  Hysteresis  for  the  1940  El  Centro 

North  South  Earthquake  Ground  Motion 


Lateral  Hysteresis  of  Corner  Node 
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Lateral  Corner  Node 


Figure  64.  30-Story,  Maxwell  Isolator,  Comer  Node  Lateral  Displacement  History  for 
the  1940  El  Centro  North  South  Earthquake  Ground  Motion 


Comer  Node  Uplift 


Figure  65.  30-Story,  Maxwell  Isolator,  Comer  Node  Uplift  Displacement  History  for  the 


1940  El  Centro  North  South  Earthquake  Ground  Motion 


Again,  Table  10  and  Table  11,  present  the  computational  results  of  the  synthesis 


process  using  MATLAB  5.2  on  a  733  Pentium  III,  Dell  Computer. 


Blocks 

Element 

2 

4 

6 

8 

Spring 

1388 

823 

783 

870 

Real  Bilinear  Element 

840 

611 

615 

782.7 

Ideal  Bilinear  Element 

847 

402 

489 

625 

Wen  Element 

819 

576 

785 

989 

Maxwell  Element 

603 

485 

576 

680 

Table  10.  30-Story,  1940  El  Centro  Earthquake  Recursive  Block-By-Block  Synthesis 

Time  Results  in  seconds 


Blocks 

Element 

2 

4 

6 

8 

Spring 

1.21E+11 

4.50E+10 

2.91  E+10 

2.34E+10 

Real  Bilinear  Element 

6.80E+10 

3.23E+10 

2.27E+10 

2.06E+10 

Ideal  Bilinear  Element 

3.15E+10 

1.97E+10 

1.74E+10 

1.60E+10 

Wen  Element 

9.13E+10 

4.50E+10 

2.90E+10 

2.33E+10 

Maxwell  Element 

3.82E+10 

2.56E+10 

2.11E+10 

1.82E+10 

Table  11.  30-Story,  1940  El  Centro  Earthquake  Block-By-Block  Synthesis  Flop  Count 

Results 

For  the  above  calculations,  the  cutoff  frequency  was  set  to  1.0  Elz,  which  retained 
a  total  of  26  modes. 

Looking  at  the  optimum  number  of  blocks  versus  time,  Table  12,  tabulates  the 
fastest  times  for  the  varying  numbers  of  isolators  and  DOFs.  The  table  also  compares 
time  synthesis  for  various  computers  and  processor  speeds.  The  synthesized  values  are 
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taken  from  a  50  second  time  duration  of  the  1940  El  Centro  earthquake  displacement 
data  for  the  nonlinear  Wen  element. 


9-Bay  4-Story 

1 5-Bay  30-Story 

8 

32 

162 

Modes  retained 

20 

20 

26 

NDOF 

120 

480 

15066 

Optimum  Number  of 
Blocks 

8 

5 

4 

10 

8 

8 

iiii 

15 

60 

1145 

Dell  Pentium  III  750 
MHz  Time  fe} 

CO 

38 

576 

Table  12.  Optimum  Times  for  Nonlinear  Wen  Element  Synthesis 

As  Table  12  indicates  that  for  fairly  large  base  isolated  structures  the  optimum 
number  of  blocks  doesn’t  deviate  much  from  smaller  models.  Although  relatively 
constant  the  timesavings  can  be  significant  from  the  original  integral  method. 
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VI.  CONCLUSIONS 


The  recursive  block-by-block  method  is  considerably  efficient  in  computing  the 
seismic  transient  dynamic  response  of  multiple  degrees  of  freedom  structures  using 
hysteretic  type  isolators.  As  Table  12  outlines,  there  are  limits  to  the  number  of  blocks 
used  to  save  computation  time.  Though  the  number  of  iterations  is  generally  reduced  per 
block,  the  recursive  iterative  process  still  requires  a  minimum  amount  of  iterations  to 
converge.  Therefore,  by  increasing  the  number  of  blocks  will  incur  added  synthesis  time 
even  though  the  number  of  floating  point  operations,  FLOPS,  is  decreased. 

Although  the  initial  reductions  in  time  computation  are  attributed  to  the  retention 
of  significant  modes  and  computing  the  synthesis  over  the  retained  (Cset)  degrees  of 
freedom,  it  is  the  recursive  block-by-block  method  that  can  further  reduce  it  up  to  an 
additional  60%.  For  large  (Cset)  models  the  major  computational  elements  are  the 
impulse  response  matrix  assembly  and  the  recursive  iterative  synthesis.  As  stated  above 
the  iterative  synthesis  will  reach  a  minimum,  however,  the  impulse  response  function 
program  can  contain  a  significant  amount  of  loops  which  can  have  a  moderate  effect  on 
time.  However,  if  the  time  duration  and  the  number  of  modes  retained  remain  constant  a 
flag  variable  can  be  entered  as  to  skip  this  computation  if  the  only  thing  that  changes  is 
isolator  properties.  This  can  significantly  reduce  the  computation  time  for  large  number 
of  isolator  based  structures. 

Since  the  program  only  iterates  on  the  retained  (Cset)  nodal  degree  of  freedoms, 
it  will  produce  the  transient  dynamic  response  history  and  the  respective  synthesized 
force  time  history.  Because  of  this,  if  two  structures  had  an  equal  number  of  base 
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isolators  and  the  same  number  of  modes  retained,  the  program’s  run  time  will  be  similar 
even  though  structure  two  had  10  times  as  many  nodes  as  structure  one  did.  To  compute 
non-isolator  coupled  degrees  of  freedom  time  history  responses  the  integral  equation 
formulation  of  the  synthesized  forces  and  the  associate  impulse  response  function  is 
employed.  This  integral  calculation  can  be  computed  directly  or  by  the  block-by-block 
method.  Either  method  produces  the  responses  quickly  and  in  a  matter  of  seconds. 

The  force  extraction  of  an  isolated  structure  using  the  trapezoidal  integration 
scheme  proved  very  accurate.  However,  the  limitations  lay  in  MATLAB’s  memory 
intensive  backslash  command,  CV,  and  the  process  only  works  on  system  that  has  only 
one  node  in  vibration. 
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VII.  RECOMMENDATIONS  FOR  FUTURE  WORK 


The  work  contained  herein  was  conducted  to  develop  and  validate  a  recursive 
block-by-block  program  that  uses  nonlinear  hysteretic  type  isolators  for  base  isolation. 
There  have  been  studies  conducted  on  the  exact  nature  of  uplift  effects  on  isolators.  The 
shared  data  from  these  studies  or  additional  laboratory  data  can  be  incorporated  into  the 
program  as  data  files,  which  MATLAB  can  easily  interpolate  from.  Using  this  data,  a 
model  can  be  designed  and  isolated  using  a  known  or  fabricated  isolator.  This  model  can 
then  be  used  for  further  validation  of  the  program.  Additionally,  optimizing  the  program 
code  can  greatly  reduce  the  computation  time  by  re-writing  the  loops  into  a  more 
efficient  manner,  specifically,  in  assembling  the  impulse  response  function.  Finally,  an 
optimization  method  can  be  employed  to  give  the  designer  an  optimal  isolator  for  a  given 
excitation  and  set  of  structural  limitations. 
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APPENDIX:  TRAPEZOIDAL  FORCE  EXTRACTION  VERSUS 


NASTRAN  AND  ODE45 

A  W36x848  50  ksi  steel  column  was  modeled  using  16  nodes  spaced  equally 
over  the  total  length  of  the  column  of  375  feet.  The  column  was  chosen  from  the 
designed  30  story  building.  The  MATLAB  code,  developed  for  the  trapezoidal  force 
extraction’s  data  was  compared  to  the  large  mass  method  used  in  NASTRAN,  [Ref.  23]. 
NASTRAN  generated  the  modal  data  as  well.  The  excitation  was  conducted  in  the  weak 
axis  of  bending  of  the  column.  A  9-DOF  system  was  used  to  validate  the  extraction 
process  against  MATLAB ’s  ode45  function. 

The  validation  is  conducted  with  a  sinusoidal  input  with  an  amplitude  of  1 
centimeter  converted  to  inches  with  a  frequency  of  .5  Hz.  In  NASTRAN,  the  time  step 
was  conducted  at  .1  seconds  for  a  duration  of  20  seconds.  Figure.  66,  depicts  the 
transient  time  history  of  the  tip  of  the  beam.  Conducting  the  Trapezoidal  force  extraction 
program  at  a  time  step  of  .01  seconds  and  a  duration  of  20  seconds  and  a  large  mass 
representing  the  ground,  yields  a  solution  that  is  very  similar  to  the  one  generated  by 
NASTRAN.  A  better  representation  of  the  two  plots  would  be  to  present  the  relative 
difference,  as  shown  in  Figure  67.  Because  the  time  response  crosses  the  time  axis, 
calculating  the  relative  errors  vice  differences  would  generate  gross  errors.  As  the  plot 
indicates  the  comparisons  are  acceptable. 
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Figure  66.  NASTRAN  Sinusoidal  Base  Excitation  Tip  Displacement  Results  of  a  375 

Feet  Steel  W36x487  Column 


Figure  67.  375  Feet  Steel  W36x487  Column ,  NASTRAN  Versus  Trapezoidal  Time 

History  Difference  Plot  of  Tip  DOF 


Using  a  9-DOF  spring  mass  system  as  presented  in  Figure  16,  MATLAB’s  ode45 
function  was  used  in  accordance  with  the  following  state  space  representation  of  the 
governing  second  ordinary  differential  equation: 

H  M  ,  ,  (62) 

j sj  =  -[MY1  [C]|x|  -  [M]-'  [. Kjx }  +  [MY]  {if0 } 

Where: 

/r=-oo— )wi  )-xo/n  (63) 

in 

Using  Figure  17,  the  trapezoidal  extraction  method  was  used  using  a  large  mass 
to  represent  the  ground  and  conducted  at  time  step  of  .01  for  10  seconds.  Figure  68, 
illustrates  the  time  history  response  of  the  first  mass.  Again,  the  ode45  plot  is  very 
similar  and  a  difference  plot  was  graphed  in  Figure  69.  From  the  figures  the  extraction 
method  is  an  acceptable  method  in  developing  a  transient  response. 
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Figure  69.  9-DOF  Spring  Mass  System,  ODE  Versus  Trapezoidal  Extraction  Absolute 

Difference  Plot  of  Base  DOF 
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